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Abstract 

Brain metabolism is controlled by complex regulation mechanisms. As part of their nature many 
complex systems show scaling behavior in their timeseries data. Corresponding scaling exponents 
can sometimes be used to characterize these systems. Here we present maps of scaling exponents 
derived from BOLD functional Magnetic Resonance Imaging (fMRI) data, and show that they 
reveal activation patterns in the human brain with very high precision. In contrast to standard 
model-based analysis we use no prior knowledge on the experimental stimulation paradigm for 
extracting activation patterns from fMRI timeseries. We demonstrate that mental activity is one- 
to-one related to large fractal exponents, or equivalently, to temporally highly correlated processes. 
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Functional magnetic resonance imaging (fMRI) based on Blood Oxygenation Level De- 
pendent (BOLD) signal- changes allows assessment of brain activity via local haemodynamic 
variations over time [I]. It provides the highest combined spatial and temporal resolutions 
presently available for non-invasive functional brain mapping. In fMRI measurements slices 
of the human brain are sampled in a repetitive manner with mm 3 sized voxels as the smallest 
information-carrying spatial element. These represent aggregated and collective information 
on the physiological status of a large number of neurons with a temporal resolution of tens of 
ms. In a typical fMRI experiment external (e.g. visual) stimuli are presented in intervals of 
several seconds (~ 0.1 Hz), causing a change in voxel-signal intensity, delayed and blurred 
by the haemodynamic response [fj, ||. Functional MRI voxel-timeseries exhibit consider- 
able noise components which adversely effect conventional fMRI analysis and therefore are 
typically considered undesired and interfering. On the other hand, by studying fMRI noise 
evidence has been found that voxel-timeseries potentially exhibit 1/f noise characteristics 
0, 0. 1/f noise - e.g. characterized by a power-law decay of the Power Spectrum Density 
(PSD) - has been discovered in a wide variety of systems which are governed by complex 
physiological regulation mechanisms such as the cardiovascular system || 0] . 

Brain metabolism is controlled by a variety of complex mechanisms that guarantee 
proper brain function including sufficient blood supply to meet neural demands. Increased 
metabolic demand is triggered by mental processes and initiates a cascade of biochemical 
reactions, followed by changes in haemodynamic parameters such as blood flow and oxy- 
genation H . In this work we argue that the use of high temporal- resolution BOLD fMRI in 
combination with a scaling analysis enables us to distinguish different physiological states 
(active and non-active) of the brain. We show that this is a robust method and can be seen 
already by using relatively simple fractal noise-measures like Hurst- and PS'D-exponents. 
We focus on the statistical scaling structure of voxel-timeseries obtained in fMRI experiments 
and visualize the amount of information present in the "noise" . 

Noise structure can be characterized by so-called fractal measures. Usually one uses more 
than one method to avoid well known systematic deficiencies inherent to either method j9]. 
Here we use two of the simplest methods available. 

Hurst Exponent: The concept of Brownian motion (random walk) can be generalized to 
fractional Brownian motion [ID| by introducing a scaling exponent H, often called Hurst 
exponent. Maybe the most transparent way to define it is by the two-point correlation 
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function of a (stochastic) data timeseries x(t), 

C(r) = ([x(t + r)-x(t)] 2 )ocT 2H , (1) 

with H being a real number, < H < 1. H serves to characterize the process in terms of 
temporal correlations: while for classical Brownian motion (H = 0.5) there are no correla- 
tions, for H > 0.5 the process is called positively correlated or persistent, i.e., if the process 
was moving upward (downward) at time t it will tend to continue to move upward (down- 
ward) at future times t' > t as well. This means that increasing (decreasing) trends in the 
past imply - on average - increasing (decreasing) trends in the future |Tl| . If H < 0.5 the 
corresponding process is called anti-correlated or anti-persistent, meaning that increasing 
(decreasing) trends in the past imply - on average - decreasing (increasing) trends. 

Power-Spectrum exponent: Alternatively, one can quantify internal noise structure in 
stochastic processes by the power-spectrum exponent S, defined by the power characterizing 
the decay of the (discrete) PSD, i.e., the squared Fourier spectrum of the process x(t), 
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oc u~ s , (2) 

with N being data size. For perfect 1-dimensional fractal processes it can be shown that H 
and S are related by S = 2H + 1. Hence, exponents H (S) different from a value of 0.5 (2) 
indicate a self affine, internal structure of the timeseries. 

Functional MRI was performed on a 3 Tesla whole-body MR scanner (MedSpec S300, 
Bruker, Ettlingen, Germany) using single-shot gradient echo-planar imaging (EPI), with a 
Tr of 200 ms and an effective Te of 28 ms. Four 4 mm thick slices parallel to the calcarine 
sulcus were measured using a field-of-view of 29 x 29 cm 2 and a matrix size of 64 x 64 pixels. 
This setup was repeated 1500 times resulting in a total measurement time of 5 minutes. 
A visual stimulation paradigm was presented at randomly distributed time points during 
scanning, separated by an interval of at least 10 seconds, guaranteeing that all stimulation 
frequencies are below 0.1 Hz. Stimulation consisted in presenting a rotating checker-board 
disk to the subjects for three seconds via video projector and mirror. 

Voxel time-course information I(x, t) was obtained from images taken at time instances t 
(x is voxel-position in 3D space). For further analysis we subtracted the temporal mean from 
I(x,t) to obtain the process I(x,t) = I(x,t) — (I(x,t)) t . This process serves as the starting 
point for all further analysis, no more pre-processing steps, neither filtering, normalization, 
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nor motion correction were applied to the data. For fractal analysis we use the cumulative 
voxel-intensity process x(x, t) = Y^t=i 

Standard fMRI analysis of the data sets was performed using linear regression, based 
on the known stimulation paradigm convolved with a model haemodynamic response func- 
tion, resulting in a model function HR(t) |T2|] . The results of this analysis are r-maps, 
i.e. correlation coefficients of the model function HR(t) with actual fMRI data for all in- 
dividual pixels at positions x. Results are thresholded to reveal the pixels activated during 
visual stimulation. In the following a t-threshold of 5 was chosen which corresponds to a 
Bonferroni-corrected p-value of 0.001, representing highly significant activation. 

Fig. H displays the averaged time-course of all activated pixels of a single, 4-slice exper- 
iment as found by standard regression analysis fll3" |. The corresponding functional r-maps 



in Fig. [|a clearly show the expected functional activation in the primary visual areas. The 
respective time-course is in excellent agreement with the stimulus applied, as indicated by 
the red bars in Fig. |l|. In Fig. |2|b and c, maps of H and S values are shown for all voxels 
within the brain slices, respectively. The color table is the same for all three parameters, 
however, adjusted to match the respective range of values. Regions with a large Hurst-value 
(red) coincide with areas of visual activation, as detected by the regression analysis with 
high statistical significance (p < 0.001). The same regions can be found in the S'-maps, 
however, to somewhat less spatial precision. Signals in frontal regions are caused by move- 
ment artifacts. We checked this by confirming that they reduce when applying standard 
realignment software. 

To demonstrate that the observed changes in the scaling parameters H and S are in fact 
properties of noise, we performed the following checks: (i) Residual analysis: We pixel- wise 
fitted our data to results of the model function HR(x, t) obtained from a linear standard 



fMRI model 12, 13 



I(x,t) = HR(x,t)+£(x,t) , (3) 

where I(x,t) is the fMRI signal and £(x,t) the residual noise at pixel position x. We 
performed the same scaling exponent analysis as before starting from the residual dataset 
£(x,t). As expected the correlation coefficients become vanishingly small and can not serve 
as a reliable measure anymore, Fig. KJa. With £(x, t) as the input, H and S are shown in Fig. 
Rib and c. The same areas of activity are recovered. This is a clear sign that £ (x, t) which is 



often (wrongly) considered "white noise" still carries relevant information. One could argue 
here that if this is the case the model is not yet good enough and has to be improved before 
discussing statistical significance of activation patterns. However, actually used models are 
based on linear approaches and can not account for the scaling effects reported here, (ii) 
Low pass filtering: By applying successive temporal smoothing of the pixel timeseries we 
find that the scaling effects eventually disappear, and obviously the structure observed here 
is a characteristics of the noise and not of the low frequency stimulation (signal), (Hi) 
Background noise is white: We find that the background noise (measurement or machine 
errors) in the signal as observed outside the head fluctuates around H = 0.5 (mean and 
standard deviation: H = 0.492 ± 0.046), as expected. 

In Fig. |] the timeseries I(x, t) for two selected voxels are shown, identified as either non- 
activated (a) or activated (b) in the standard fMRI analysis. Periods of visual stimulation are 
again marked with red bars. Clearly, the voxel in the visual region follows the stimulus (much 
like the entire cluster; compare Fig. [I]). In Fig. £|c and d we show the PSD of x(x, t). Voxel 
timeseries show typical 1/f noise characteristics and the power spectral density exponent S 
can be estimated safely by fits to the PiS-D-slope in a log- log plot. Fits are given by red 
lines. In order to focus exclusively on the noise structure and to exclude any influence of 
the BOLD response itself, the frequency range corresponding to the stimulus presentation 
(between 0.04 - 0.1 Hz) was excluded from the fit-range. Frequency fluctuations larger than 
1 Hz were also not included as those frequencies are above the fractal onset frequency for 
the given data size ||. In Fig. [|e and f we plot C(r) from Eq. [l[ and extract the Hurst 
exponent H as half of the slope of the linear fit. Again, the same effective fit-range (red 
lines) was used. 

A priori, noise components above the characteristic stimulus frequencies could be thought 
to appear as measurement noise and undesirable artifacts. In this work we have examined 
these noise components from a scaling point of view. We found that in non-active cortical 
regions, voxel-activity is well described by classical Brownian motion (random walk model, 
H ~ 0.5 and S ~ 2), while noise components from voxels of active brain regions - besides 
obviously following the stimulus - change their character towards highly correlated stochastic 
processes, i.e., correlated fractional Brownian noise. We have also observed that the latter 
regions are detected by H and S at very high sensitivity, i.e., regions of large H and closely 
correlate with large r-values from model-based regression analysis. In areas outside the 
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subject's head, H is fluctuating around 0.5 ensuring that no internal noise structure is arising 
or due to measurement or machine noise. These findings clearly demonstrate persistent 
fractal dominance of noise in voxel-timeseries in regions of neural activation, while in other 
areas uncorrelated noise prevails, consistent with recent work [TJ], [T5], [TO] . S and H seem to 
be related to some extent, indicating the voxel-processes to be close to perfect fractals. 



We find that even residual timeseries, as e.g. obtained by SPM software [13 , show 
practically the same scaling information content. This finding has far reaching consequences 
for statistical analyses of fMRI data. Standard statistics is based on models which leave 
white noise as a residuum. Based on white noise, null hypotheses can be tested and sensible 
^values be quoted. For correlated or scaling processes, however, it is known that standard 
statistics is not valid, since residua are by definition correlated (e.g. via scaling relations). 
As a consequence other methods based on more complicated statistics have to be used to 
arrive at "p-values" which capture actual statistical significance. Maybe Tsallis entropy 



17] and its associated non-extensive statistics offer a natural starting point for consistently 



solving this problem. 

In conclusion we find that the statistical nature of voxel-timeseries depends on the under- 
lying level of neural activity. In activated regions the noise processes are strongly temporally 
correlated. The practical potential of this study is that scaling analysis is entirely model- 
independent: we do not make any assumptions about the nature of the voxel-processes, but 
rather try to clarify it. Moreover, no information on the stimulus paradigm is required. It 
is intriguing that the information carried by fMRI noise alone is sufficient to at least com- 
pete with the gold-standard of fMRI analysis for detecting statistically significant regions of 
neural activation. Scaling measures - as a complementary source of information - may there- 
fore improve specificity in fMRI data analysis. We believe that the use of fractal parameters 
in fMRI also represents a significant potential for understanding the fundamental processes 
involved in neuro-vascular coupling, by identifying "noise components" as valuable measure- 
ments representing detailed high(er) frequency information, rather than treating "noise" as 
a nuisance. Our results can be seen in the context of recent trends which focus on the 



practical role of noise |]18| , and that one should not restrict oneself to the obvious, but also 
hunt for the information hidden in the "noise" . 
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FIG. 1: Average time-course of all activated pixels (blue line). Overlayed (red dashed line) is the 
functional paradigm convolved with the haemodynamic response (model function HR(t)). The 
stimulation periods are indicated by red bars at the bottom. 

FIG. 2: Computed "activation" maps of four brain slices using the standard linear regression 
analysis. Correlation coefficients (r-map) of the actual data and the functional paradigm convolved 
with the haemodynamic response function are shown (a). Hurst exponent analysis for the pixel 
timeseries (b); Power spectral exponent analysis for the pixel timeseries (c). 

FIG. 3: Same as previous figure, but starting from the residual dataset £(x,t). Clearly the corre- 
lation map (a) does not show activation patterns anymore, notice the change in scale compared to 
the previous figure. H and S maps, (b) and (c), show the same activated regions and about the 
same numerical values as in the original data. 

FIG. 4: Timeseries of a non-active (a) and an active (b) voxel. The power spectral exponent S for 
these processes is given by the slope of the PSD shown in (c) and (d). Clear 1/f behavior in the 
PSD is observable. The Hurst exponent H is defined as half the slope of the function C. 



8 




FIG. 1 

9 
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